1. Field of the Invention
The invention relates to a method of imaging an interior of a turbid medium, including consecutive irradiation of the turbid medium with light and measurement of the light propagated along a plurality of light paths through the turbid medium, and determination of an image of the interior of the turbid medium from the intensities measured. The invention also relates to a device for carrying out such a method.
2. Description of the Related Art
In the context of the present Patent Application light has to be understood as electromagnetic radiation having a wavelength in the visible or infra-red range approximately between 400 and 1000 nm. Furthermore, a turbid medium has to be understood as a volume of a highly scattering substance, for example, an Intralipid solution or biological tissue.
Such a method is known from the article "Back-projection Reconstruction of Cylindrical inhomogeneities from Frequency Domain Optical Measurements in Turbid Media", S. A. Walker et al, OSA Topics on Advances in Optical Imaging and Photon Migration, 1996. In the cited article, the known method is used for the imaging of a turbid medium having optical properties similar to those of biological tissue. In medical diagnostics the known method and device could be used for imaging the internal structure of breast tissue of a human female. For example, a tumor can be localized in such an image of the internal structure of breast tissue. The known method employes a back-projection method similar to X-ray computer tomography to obtain an image of the region of interest inside the turbid medium from physical parameters of the measured intensities. In said back-projection method, the physical parameters are calculated from frequency domain data. Furthermore, a value of the optical absorption coefficient .mu..sub..alpha. and a value of the reduced scattering coefficient .mu.'.sub.s are used as reconstruction parameters. Said values are calculated for each measurement and, in a single step, are projected back along the region between source and detector while employing an appropriate weighting function. An image is formed by averaging data from multiple source detector scans taken at multiple angles. A weighting function can also be chosen in order to maximize resolution which weighting function consists of an evenly weighted straight line between source and detector multiplied by a filtering function. That filtering function accounts for the size of the sampling interval. A drawback of the known method is that the known method is not capable of reconstructing inhomogeneities beyond a source detector line.
It is an object of the invention to provide an image reconstruction algorithm for imaging inhomogeneities also beyond the source detector line. It is a further object to provide a real-time imaging technique for use in medical imaging processes, for example mammography. Therefore, the method according to the invention is characterized in that it includes the following steps:
c) determination for each measurement of a difference between an expected photon fluence and a measured photon fluence from the intensity measured, PA1 d) determination of possible strengths assignable to each pixel of the image from combinations of weighting functions and the differences determined, PA1 e) determination of a distribution function of the possible strengths assignable to each pixel, PA1 f) determination of an image from the distribution function determined. PA1 determining intensities due to a perturbation function for the position of a strong object, PA1 subtracting the intensities determined from the intensities measured, PA1 repeating the steps c) to f). PA1 a) a light source for irradiating the turbid medium, PA1 b) a plurality of detectors for measuring light intensities of light propagated through the turbid medium at a plurality of different positions, PA1 c) means for irradiating the turbid medium from a plurality of different positions connected to the light source, PA1 d) means for selecting a detector from the plurality of detectors, PA1 e) a control unit for generating control signals for the means for controlling the means for irradiating the turbid medium and the means for selecting detectors, and PA1 f) a reconstruction unit for reconstructing an image from the measured intensities, the reconstruction unit also being arranged to carry out the following steps: PA1 1) consecutive irradiation of the turbid medium with light and detection of the light propagated along different light paths through the turbid medium, PA1 2) determination of an image of the interior of the turbid medium from the measured intensities, PA1 3) determination, for each measurement, of a difference between an expected photon fluence and a measured photon fluence from the intensity measured, PA1 4) determination of possible strengths assignable to each pixel of the image from combinations of weighting functions and the differences determined, PA1 5) determination of a distribution function of the possible strengths assignable to each pixel, PA1 6) determination of an image from the distribution function determined.
In the present Patent Application the expected photon fluence is defined as ##EQU1## in which S.sub.O represents the source intensity, K=(.mu..sub..alpha. c/D').sup.1/2 .about.(3.mu..sub..alpha. .mu.'.sub.s).sup.1/2, D'=cD=c/3(.mu.'.sub.s +.mu..sub..alpha.)! is the photon diffusion coefficient (in m.sup.2 /s) and c is the speed of light in the turbid medium (in m/s). Furthermore, the strength q of an inhomogeneity at a position r can be regarded as the effect on the measured intensity of the light an is given by the formula ##EQU2## wherein ##EQU3## and r.sub.l,r.sub.s,r.sub.d represent position vectors of a point in the object space, a light source position, and a detector position, respectively and .phi.(r.sub.d) is the measured photon density.
Furthermore, a weighting function is defined as a function giving a dependency of the strength of an object of a position r, which strength introduces a normalized intensity change at a detector position r.sub.d in light from a source at position r.sub.s. The weighting function W is defined as the inverse function of the perturbation function P, so ##EQU4## The distribution function is defined as a function giving a correlation between the possible strengths assignable to each pixel of the reconstructed image for the differences determined.
The effect of said steps is that an image is reconstructed containing an non-blurred image of a strong object. The invention is based on the insight that the possible strengths assignable to a pixel can be derived from a combination of the weighting function and the differences between the measured photon fluence and an expected photon fluence for each light source position and each detector position and that furthermore the possible strengths determined for different light source positions and detector positions must correlate in order to assign a value of the strength to that pixel. Said correlation can be determined by employing the distribution function. Consequently, determination of an image from the distribution function can be performed by employing statistics from the distribution function. Furthermore, said method enables three-dimensional imaging because the reconstruction of an image is not limited to inhomogeneities in a source detector plane only, but because of the use of the weighting function, also accounts for inhomogenities beyond the source detector plane.
A further advantage of the method is that fewer measurements of different source detector positions are suffice as compared to the known back-projection method.
An embodiment of the method according to the invention is characterized in that the method also includes the following steps:
The perturbation function P(r.sub.l,r.sub.s,r.sub.d) is defined as the inverse function of the weighting function. The intensity at a detector position is then given by I=A.sub.O .PHI.(r.sub.d)=A.sub.O .PHI..sub.O (r.sub.d)1-qP(r.sub.l,r.sub.s,r.sub.d)!, wherein A.sub.O (in cm.sup.2) a collection area of the input window f the detector, q represents the strength of the perturbation and .PHI..sub.O(r.sub.d) is defined as ##EQU5## in which S.sub.O represents the initial source intensity and ##EQU6## represents the normalized photon diffusion coefficient. A strong object is an indication for an inhomogeneity having a large value of q. A small object is an indication for an inhomogeneity with a small value of q.
Subtracting the effect of the strong object from the intensities measured and repeating said steps of the image reconstruction process produces an image in which smal objects are deblurred. Furthermore, the image can be refined, by a small number of these iterations, so that in the resultant image more smal objects become visible.
An embodiment of the method according to the invention is characterized in that the distribution function comprises a histogram for each pixel, in which histograms each bin represents a value of the possible strength and the frequency of the bin represents the occurrence of the possible strengths determined. An image of the interior of the turbid medium is then formed by using the statistical information of said histogram. Generally speaking, a high frequency indicates an inhomogeneity for that pixel in the turbid medium. Statistical parameters that may used for forming the image are, for example, maximum value or peak value, mean, variance and standard deviation.
An embodiment of the method according to the invention is characterized in that for each pixel the value of the bin belonging to maximum frequency of the histogram of the pixel is assigned to the pixel. The effect of this step is that the possible strength of the pixel on which the different possible assignable strengths agree or correlate is assigned to the pixel.
An embodiment of the method according to the invention is characterized in that for each pixel the pixel value assigned is multiplied by a correction factor ##EQU7## wherein .sigma..sub.k is the standard deviation of the histogram of the pixel and a is a constant determined by the desired bit accuracy of the image as ##EQU8## wherein b represents the bit accuracy and q.sub.max is the maximum value of the weighting function. Bit accuracy is defined as the number of binary digits representing the maximum possible values of a pixel in the image. The effect of this step is that the strength of the pixel is corrected for the standard deviation of the histogram. For a high standard deviation .sigma..sub.k only a small value of the strength will be assigned to the pixel. For a low standard deviation .sigma..sub.k the histogram is peaked around a central value and that central value will be assigned to the pixel.
The invention also relates to a device, which is provided with
characterized in that the reconstruction unit is further arranged to carry out the following sub-steps for determining the image of the object:
The above and other, more detailed aspects of the invention will be explained hereinafter by way of example with reference to the drawing.